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Fast method for quantum mechanical molecular dynamics 
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As the processing power available for scientific computing grows, first principles Born- 
Oppenheimer molecular dynamics simulations are becoming increasingly popular for the study of a 
wide range of problems in materials science, chemistry and biology. Nevertheless, the computational 
cost of Born-Oppenheimer molecular dynamics still remains prohibitively large for many potential 
applications. Here we show how to avoid a major computational bottleneck: the self-consistent-field 
optimization prior to the force calculations. The optimization-free quantum mechanical molecular 
dynamics method gives trajectories that are almost indistinguishable from an "exact" microcanoni- 
cal Born-Oppenheimer molecular dynamics simulation even when low pre-factor linear scaling sparse 
matrix algebra is used. Our findings show that the computational gap between classical and quan- 
tum mechanical molecular dynamics simulations can be significantly reduced. 



I. INTRODUCTION 

The past three decades have witnessed a dramatic in- 
crease in the use of the molecular dynamics simulation 
method [U |2] ■ While it is unquestionably a powerful and 
widely used tool, its ability to calculate physical prop- 
erties is limited by the quality and the computational 
complexity of the interatomic potentials. Among com- 
putationally tractable models, the most accurate are ex- 
plicitly quantum mechanical with interatomic forces cal- 
culated on-the-fly using a nuclear potential energy sur- 
face that is determined by the electronic ground state 
within the Born-Oppenheimer approximation [2H1]- In 
Hartree-Fock [5J [6] or density functional theory [TlUO). 
the electronic ground-state density is given through a 
self-consistent-field (SCF) optimization procedure, which 
involves iterative mixed solutions of the single-particle 
eigenvalue equations and accounts for details in the 
charge distribution. Since the interatomic forces are sen- 
sitive to the electrostatic potential [11], molecular dy- 
namics simulations are often of poor quality without a 
high degree of self-consistent-field convergence. This is 
unfortunate since the iterative self-consistent-field proce- 
dure is computationally expensive and in practice always 
approximate. 

Recently there have been efforts to reduce the com- 
putational cost of the self-consistent-field optimization 
without causing any significant deviation from "exact" 
Born-Oppenheimer molecular dynamics simulations [121 - 
03] . In this article we go one step further, and in analogy 
to time-dependent techniques such as Ehrenfest molecu- 
lar dynamics [T5HT7] or the Car-Parrinello method [2jH8l- 
123] . we show how the electronic ground state optimiza- 
tion can be circumvented fully without any noticeable 
reduction in accuracy in comparison to "exact" Born- 
Oppenheimer molecular dynamics. 

Our optimization-free dynamics is based on a re- 
formulation of extended Lagrangian Born-Oppenheimer 



'Corresponding Author Email: amn@lanl.gov 



molecular dynamics |35j in the limit of vanishing self- 
consistent-field optimization. The method is presented 
within a general free energy formulation that is valid also 
at finite electronic temperatures and should be applica- 
ble to a broad class of materials. In addition to the re- 
moval of the costly self-consistent-field optimization we 
also demonstrate compatibility with low pre-factor lin- 
ear scaling electronic structure theory [24lf26] . The com- 
bined scheme provides a very efficient, energy conserving, 
low-complexity method for performing accurate quantum 
molecular dynamics simulations. 



II. FAST QUANTUM MECHANICAL 
MOLECULAR DYNAMICS 



A. Born-Oppenheimer molecular dynamics 

Born-Oppenheimer molecular dynamics based on den- 
sity functional theory can be described by the Lagrangian 



£ BO (R, R) = lYl M ^t ~ U l R -> /"]> 



(1) 



where the potential energy, 



U[K;p]=2j2 e i 



p{v)p(r>) 



dr'dr 



(2) 



V xc [p}p(r)dr + E xc [p] + £ ZZ [R], 



is calculated at the self-consistent electronic ground state 
density, p{r), for the nuclear configuration R = {Rk} 
[in] . Here, £j are the (doubly) occupied eigenvalues of 
the effective single-particle Kohn-Sham Hamiltonian, 



H[ P ] 



-V 2 + K(R,r) 



P(r') 



■dr' + V xc [p], (3) 



where V^ c [p] is the exchange correlation potential, 
Vn(R, r) the external (nuclear) potential, and —5V 2 the 
kinetic energy operator. -E ZZ [R] is the electrostatic ion- 
ion repulsion and E xc [p] the exchange correlation energy. 



If the electron density deviates from the ground state 
density p by some small amount Sp, the error in the po- 
tential energy is essentially of the order Sp 2 , depend- 
ing on the particular formulation used for calculating 
U[R;p + Sp] [331 |3S S3- However, since the Hellmann- 
Feynman theorem is valid only at the ground state den- 
sity, we do not have a simple expression for the forces 
that avoids calculating derivatives of the electronic den- 
sity, d(p + Sp)/dRk- In practical calculations, the ac- 
curacy of the potential energy can therefore not be ex- 
pected to hold also for the forces and a high degree of 
self-consistent-field convergence is therefore typically re- 
quired. 



B. Extended Lagrangian molecular dynamics 

Here we outline how we can circumvent the self- 
consistent-field procedure in Born-Oppenheimer molec- 
ular dynamics. Instead of recalculating the ground 
state density before each force evaluation with an it- 
erative optimization procedure, the idea here is to use 
an auxiliary density n(r), as in extended Lagrangian 
Born-Oppenheimer molecular dynamics J35H39) . which 
evolves through a harmonic oscillator centered around 
the ground state density p(r). Based on a general 
free energy formulation of extended Lagrangian Born- 
Oppenheimer molecular dynamics [39 in the limit of van- 
ishing self-consistent-field optimization, we define the ex- 
tended Lagrangian: 



£(R,R,n,n) = lYM k Rl 



2^ 

k 

-p / h(r) 2 dr- -puj 2 



-U[R;n]+T e S[R;n] 



(p(r) — n(r)) dr. 



While the potential and entropy terms, IA and S, are 
well defined at the ground state density [3], i.e. when 
n = p, there are several different options when n deviates 
from p, e.g. the Harris- Foulkes functional [33j [2H 147] . 
In a more general case, the potential energy and en- 
tropy term may therefore also be determined by n(r) 
implicitly through an additional function cr[n(r)], i.e. 
U[R;n] = U[R;n,a[n]] andS[R;n] = <S[R; n,cr[n]]. Here 
<r[n(r)] is a temperature dependent density given from 
the diagonal part of the real-space representation of the 
(doubly occupied) density matrix, which is given through 
a Fermi-operator expansion 9 of the effective single- 
particle Hamiltonian, i?[n], i.e. 



a(r) = a[n(r)] = 2 ( e «H[n]-w>/) + ]Y 



(5) 



At zero electronic temperature the Fermi-operator ex- 
pansion corresponds to a step function with the step 
formed at the chemical potential, pq. In our Lagrangian 
above, p and u are fictitious mass and frequency pa- 
rameters of the harmonic oscillator and /3 is the inverse 



electronic temperature, i. e. = l/{k-B,Te)- The pur- 
pose of the entropy- like term S[R; n] is here to make the 
derived forces of our dynamics variationally correct for 
a given entropy- independent density, n(r), at any elec- 
tronic temperature. This approach is different from the 
regular formulation where the density is determined by 
the entropy through the minimization of the electronic 
free energy functional [pj H0HI2"] . 



1. Equations of motion 

The molecular trajectories corresponding to the ex- 
tended free energy Lagrangian C in Eq. Q are deter- 
mined by the Euler-Lagrange equations of motion, 



oR k 


, rr, dS[R;n] 

n C 9R k 


puj 2 d f 
- 2 dRkJ {p{r) - 


- n(r)) dv 


n 



(6) 



and 



ph{v) = puj 2 (p(r) — n(r)) 



SU[R-/, 



T, 



SS\Ry, 



R 



(7) 



R 



where the partial derivatives are taken with respect to 
constant density, n, or coordinates, R. The limit p — > 
gives us the equations of motion of our extended La- 
grangian dynamics, 



M k R k 



dU[R;n] 



dR k 



T, 



dS[R; n] 



dR k 



n(r) = uj 2 (p(r) -ra(r)), 
where we have defined S[R; n] such that 



5U[R;; 



= T, 



<55 [R; 



R 



(8) 



(9) 



(10) 



R 



As is shown in the Appendix (Sec. VII A), the corre- 



sponding property for dS/dR k is also of importance for 
the calculation of the Pulay force in Eq. |8]). Notice 
that these equations still require a full self-consistent 
field optimization, since the auxiliary density n(r) evolves 
around the ground state density p(r). 

Since the nuclear degrees of freedom do not depend on 
the mass parameter p in Eqs. ^ and ^, the total free 
energy, 

k 

is a constant of motion in the limit of vanishing p. More- 
over, if Etot is close to the exact ground state free energy 



for approximate densities n(r), we can also expect that 
the forces of the extended Lagrangian dynamics should 
be accurate. 

The forces in Eq. (ISj) are calculated at the approximate, 
unrelaxed, density n(r) using a Hcllmann-Feynman-like 
expression, where the partial derivatives are taken with 
respect to a constant density n(r). This is possible only 
because n(r) appears as an independent dynamical vari- 
able. In general, as mentioned above, this can not be 
assumed, since the Hellmann-Feynman force expression 
is formally applicable only at the ground density. A more 
detailed derivation of explicit force expressions is given 
in the Appendix. 



2. Entropy contribution 

Depending on the particular functional form chosen for 
the potential energy term, U(R; n), we may not have ac- 
cess to a simple explicit expression of <S[R; n] that fulfills 
Eq. (10 1. In this case an approximate entropy term has 
to be used. This has no effect on the dynamics in Eqs. $h 
and Q, since the forces remain exact by definition. An 



approximation of the entropy term therefore only affects 
the estimate of the constant of motion, -Etoti in Eq. (111. 



We have found that the regular expression for the elec- 
tronic entropy [S], 



S[R; n] = -2fc B V {ft ln(/«) + (1 - /,) ln(l - /,)} , 



(12) 
which formally is defined only at the ground state den- 
sity, i.e. when n = p, typically provides a highly accurate 
approximation also for approximate densities as will be 
illustrated in the examples below. Here /, are the occu- 
pation numbers of the states, i.e. the eigenvalues of the 
density matrix in Eq. (T5J) . These are determined by the 
Fermi-Dirac distribution of the single-particle eigenvalues 
E; of the Hamiltonian H\n], i.e. 



f, 



= e^ 



1 



(13) 



By comparing the calculation of -Etot in Eq. (11) using 
the approximate entropy term, <S[R;n], in Eq. ( |12[ ) to 
"exact" , fully optimized, Born-Oppenheimer molecular 
simulations, we can estimate the accuracy of our dynam- 
ics. 



C. Fast quantum mechanical molecular dynamics 

As in extended Lagrangian Born-Oppenheimer molec- 
ular dynamics, the irreversibility of regular Born- 
Oppenheimer molecular dynamics that is caused by the 
self-consistent-field optimization, can be avoided, since 
the density n(r) can be integrated using a reversible ge- 
ometric integration algorithm [351 [Ml 1331 S3, e.g. the 
Verlet algorithm as in Eq. ( 16 1 below. This prevents the 



unphysical drift in the energy and phase space of regular 
Born-Oppenheimer molecular dynamics |12H14j and our 
dynamics will therefore exhibit long-term stability of the 
free energy E tot in Eq. (fill. 



A main problem so far is that we still need to calcu- 
late the self-consistent ground state density p(r) in the 
integration of n(r) in Eq. pj. Fortunately, various geo- 
metric integrations of the auxiliary density n(r) in Eq. ([9]) 
are stable also for approximate ground state density esti- 
mates of p(r), as long as the approximation of p(r) is at 
least infinitesimally closer to the exact ground state com- 
pared to n(r). Using an integration time step of St, this 
stability holds if the value of the dimensionless variable 
K = St 2 uj 2 is chosen to be appropriately small [36, 45, 46J. 



For energy functionals that are convex in the vicinity of 
the ground state density we may therefore replace p(r) in 
Eq. (|9| by a linear combination (1 — c)n + ca 43! , which 
gives us the approximate equations of motion 



M k R k = - 



du[n v , 



dR k 



■T„ 



dS[Ry, 



OR, 



and 



n(r) = oj 2 (cr(r) — n(r)), 



(14) 



(15) 



where the constant to 2 has been rescaled by c. The Verlet 



integration of Eq. (151, including a weak dissipation to 



avoid an accumulation of numerical noise [361 137] , 



n t+St = 2n t - rit-st + St w (a t - n t ) + 



K 

«y^c fc n t _ fc g t , 

fc=0 



(16) 

is therefore stable if a sufficiently small positive value 
of k — St 2 oj 2 is chosen [33] ■ Thus, without any self- 
consistent-field optimization of p(r), the previously opti- 
mized values of k in Ref. [3H1SS1H5] should be rescaled by 
a positive factor < 1. Certain ill behaved (non-convex) 
functionals with self-consistent-field instabilities [13] can 
not be treated in this framework. 

The proposed molecular dynamics as given by Eqs. 
( 14 ) and ( 15 ) is the central result of this paper. The 



equations of motion do not involve any ground state 
self-consistent-field optimization prior to the force eval- 
uations and only one single diagonalization or density 
matrix construction is required in each time step. The 
frequency w of the electronic density is well separated 
from the nuclear vibrational oscillations. Using a value 
of Stui = t/k « 1 and an integration time step St, which is 
~ 1/15 of the period of the nuclear motion, the frequen- 
cies differ by a factor of 5. As will be demonstrated in the 
examples below, the scheme is also fully compatible with 
linear scaling electronic structure theory [531 122] ■ This 
compatibility is crucial in order to simulate large sys- 
tems. The removal of the costly ground state optimiza- 
tion, in combination with low-complexity linear scaling 
solvers, provide a computationally fast quantum mechan- 
ical molecular dynamics (Fast-QMMD) that can match 
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FIG. 1: Total energy fluctuations, Eq. (Ill, using "exact" (4 



SCF/step) Born-Oppenheimer molecular dynamics (BOMD), 
and the fast quantum mechanical molecular dynamics, Eqs. 
(14} and (THJ, (Fast-QMMD), with (r > 0) or without (r = 0) 
thresholding applied in the low pre-factor linear scaling solver 
|26| . The simulations were performed with the molecular dy- 
namics program latte using self-consistent-charge density 
functional based tight-binding theory in an orthogonal for- 
mulation at T e — 0, i.e. as in Eqs. (571, (581 and (60 1. 



FIG. 2: Panel a) shows the x-plane phase space trajectory of 
a single carbon atom (C) based on an "exact" (4 SCF/step) 
Born-Oppenheimer molecular dynamics (BOMD, dashed line) 
and the fast quantum mechanical molecular dynamics (Fast- 
QMMD, solid line). Panel b) shows the fluctuations in the 
net auxiliary charge rii{t) and ground state charge qi(t) for 
the same carbon atom (i—G). The numerical threshold r 
is applied in the linear scaling solver [26]. The simulations 
were performed with the program latte using self-consistent- 
charge tight-binding theory in an orthogonal formulation at 
zero electronic temperature, i.e. as in Eqs. (571, (58 1 and (60 1. 



the fidelity and accuracy of regular Born-Oppenheimer 
molecular dynamics. 

There are several alternative approaches to derive or 
motivate the equations of motion of the fast quantum 
mechanical molecular dynamics, Eqs. (14) and (15), and 



details of the dynamics may vary depending on the choice 
of the functional form of U(R; n). However, the particu- 
lar derivation presented here is the most transparent and 
general approach that we have found so far. 

The equations of motion are given in terms of the elec- 
tron density, but they should be generally applicable to 
a large class of methods, such as Hartree-Fock theory, 
which is analyzed in the Appendix (Sec. VII A), or plane 
wave pseudo-potential methods |37) . Here we will demon- 
strate our fast quantum mechanical molecular dynam- 
ics scheme using self-consistent-charge density functional 
tight-binding theory 47-50 , as implemented in the elec- 
tronic structure code latte [IJT] , either with an orthogo- 
nal or a non-orthogonal representation and both at zero 
and at finite electronic temperatures. With this method 
we can easily reach the time and length scales neces- 
sary to establish long-term energy conservation and lin- 
ear scaling of the computational cost. Details of the com- 
putational method and our particular choices of li(R;n) 
are given in the Appendix. 



TABLE I: Wall clock timings of the fast quantum mechanical 
molecular dynamics (Fast-QMMD) simulations in compari- 
son to Born-Oppenheimer molecular dynamics (BOMD) (4 
SCF/step), without (r = 0) and with (r > 0) a low pre- 
factor linear scaling solver for the density matrix [26] with 
threshold tolerance r. The program (latte in its orthogonal 
formulation at T e = 0) was executed on a single core of a 2.66 
GHz Quad-Core Intel Xeon processor. 



Polyethene chain C100H202 


Efficiency 


BOMD (r = 0) 
Fast-QMMD (r = 0) 
Fast-QMMD (r = 10~ 5 ) 


7.5 s/step 
1.5 s/step 
0.61 s/step 


Liquid Methane (CH4)ioo 


Efficiency 


BOMD (r = 0) 
Fast-QMMD (r = 0) 
Fast-QMMD (r = 10~ 5 ) 


12.5 s/step 
2.5 s/step 
0.35 s/step 



III. EXAMPLES 

A. Orthogonal representation 

Figure [l] shows the fluctuations in the total energy (ki- 
netic plus potential) using the fast quantum mechani- 
cal molecular dynamics, Eqs. (14) and (15), as imple- 
mented in Eqs. ( [57| , (58) and (60), and an "exact" 
Born-Oppenheimer molecular dynamics |35j . for liquid 
methane (density = 0.422 g/cm 3 ) at room temperature. 



The calculations were performed with the latte molec- 
ular dynamics program using periodic boundary condi- 
tions and an integration time step of St = 0.5 fs. Since 
the molecular system is chaotic, any infinitesimally small 
deviation will eventually lead to a divergence between 
different simulations. However, even after hundreds of 
time steps and over 300 fs of simulation time the total 
energy curves are virtually on top of each other as is seen 
in the inset. The same remarkable agreement is seen in 
Fig. [2] which shows the projected phase space of an indi- 
vidual carbon atom and the fluctuations of its net charge. 
In this case the C atom was displaced compared to the 
simulation in Fig. [I] to further enhance the charge fluc- 
tuations. 
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B. Linear scaling 



The fast quantum mechanical molecular dynamics 
scheme is also stable in combination with approximate 
linear scaling sparse matrix algebra [2H US] . Using the 
recursive second order spectral projection method for the 
construction of the density matrix |26j with a numerical 
threshold, r = 10 -5 , below which all elements are set to 
zero after each individual projection, we notice excellent 
accuracy and stability in Fig. IT] without any systematic 
drift in the total energy. 

Despite their high efficiency and low computational 
pre-factor compared to alternative linear scaling elec- 
tronic structure methods [52] . it has been argued that 
recursive purification algorithms are non-variational and 
therefore incompatible with forces of a conservative sys- 
tem |25j . which is necessary for long-term energy con- 
servation. As is evident from Figs. [I] and [2] this is 
not a problem. The graphs are practically indistinguish- 
able from "exact" Born-Oppcnheimcr molecular dynam- 
ics, without any signs of a systematic drift in the to- 
tal energy. The corresponding linear scaling compati- 
bility with microcanonical simulations was recently also 
demonstrated for sclf-consistent-field-optimized extended 
Lagrangian Born-Oppenheimer molecular dynamics |27j . 

The gain in speed using the fast quantum mechani- 
cal molecular dynamics scheme in comparison to Born- 
Oppenheimer molecular dynamics is illustrated by the 
wall-clock timings shown in Tab. Ill 



FIG. 3: Total energy fluctuations, Eq. (Ill, using "exact" (4 



SCF/step) Born-Oppenheimer molecular dynamics (BOMD) , 
and the f ast quantum mechanical molecular dynamics, Eqs. 
(fli} and [l5], (Fast-QMMD), with (r = 1CT 5 ) or without 



(r = 0) thresholding applied in the low pre-factor linear scal- 
ing solver [26]. The simulations were performed with the 
molecular dynamics program latte in the non-orthogonal 
formulation at k^T B = eV, i.e. as implemented in Eqs. (50 1, 



(511 and ( 55 1 with the entropy term approximated by S — 0. 



to zero, fceTe = eV, and for the examples in Fig. [4] 
fceTe = 2 eV. In the first time step the initial nuclear 
temperature, T; n i t , was set to 300 K using a Gaussian 
distribution of the velocities. Despite the approximation 



of p in Eq. ( 15 ) and the approximate estimate of the en- 



tropy contribution to the free energy there is virtually no 
difference seen between the fast quantum mechanical and 
the Born-Oppenheimer molecular dynamics simulations. 
As in the orthogonal case, the non-orthogonal formula- 
tion of our fast quantum mechanical molecular dynamics 
is fully compatible with linear scaling complexity in the 
construction of the density matrix at T e = K. In Fig. 
[3] the reduced complexity simulation shows no significant 
deviation from "exact" Born-Oppenheimer molecular dy- 
namics. At finite electronic temperatures, the linear scal- 
ing construction of the Fermi operator [251 H2] has not 
yet been implemented. 



C. Non-orthogonal representation 

For non-orthogonal representations at finite electronic 
temperatures, a Pulay force term and a finite approxi- 
mate entropy contribution to the total free energy have 
to be included. Figures [3] and [4] illustrate the total energy 
fluctuations for the fast quantum mechanical molecular 
dynamics simulations of a hydrocarbon chain as imple- 
mented in LATTE using Eqs. (501, (51) and ( [55J ), with 
the approximate entropy term in Eq. (56). The elec- 
tronic temperature of the examples in Figure [3] is set 



D. Long-term stability and conservation of the 
total energy 

To assess the long-term energy conservation and the 
stability we use a test system comprised of 16 molecules 
of isocyanic acid, HNCO, at a density of 1.14 g cm~ 3 . 
The system was first thermalized to a temperature of 
300 K over a simulation time of 12.5 ps by the rescal- 
ing of the nuclear velocities. The simulations used an 
integration time step, St, of 0.25 ps. The simulations 
were performed using self-consistent tight-binding the- 
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FIG. 4: Total free energy fluctuations, Eq. (Ill, using "ex 



act" (4 SCF/step) Born-Oppenheimer molecular dynamics 
(BOMD), and the fast quantum mechanical molecular dy- 
namics, Eqs. fil l and |l5| , (Fast-QMMD). The simulations 
were performed with the molecular dynamics program latte 
using the non-orthogonal formulation at an electronic temper- 



ature of SibIe = 0.5 eV, i.e. as implemented in Eqs. ( 50 1 



and ( 55 1 with the entropy term approximated by Eq. ( 56 



51) 



ory [47H50J with a non-orthogonal basis as implemented 
in latte, using Eqs. ( 50 ) , ( 51 ) and ( 55 1 with the entropy 



term approximated by Eq. ( 56 ) . 

Fast quantum mechanical molecular dynamics and 
"exact" Born-Oppenheimer molecular dynamics simula- 
tions with 4 self-consistent field cycles per time step were 
performed over 250,000 time steps (62.5 ps) with T e = 
K and k&T e = 0.5 eV. The latter temperature is small 
with respect to the HOMO-LUMO gap of HNCO, which 
is about 6.0 eV, yet the entropy term, Eq. (12 1 or Eq. 



(561, contributes about 0.19 eV to the total energy ow- 



ing to the partial occupation of states in the vicinity of 
the chemical potential. Trajectories computed at T e = 
K with "exact" Born-Oppenheimer molecular dynamics 
and the fast quantum mechanical molecular dynamics 
method without (r = 0) and with (r = 10~ 5 ) linear scal- 
ing constructions of the density matrix are presented in 
Fig. [5j The standard deviation of the fluctuations of the 
total energy about its mean and an estimate of the level 
of the systematic drift of the total energies are presented 
in Table [TTJ These data show that the fast quantum me- 
chanical molecular dynamics simulations yield trajecto- 
ries that are effectively indistinguishable from the "ex- 
act" Born-Oppenheimer trajectories. Moreover, as was 
seen above, the fast quantum mechanical molecular dy- 
namics scheme appears to be fully compatible with linear 
scaling construction of the density matrix and the result- 
ing approximate forces, since this trajectory differs from 
the "exact" Born-Oppenheimer molecular dynamics tra- 
jectory only by a small-amplitude random-walk of the 
total energy about its mean [37] ■ The systematic drift 



FIG. 5: Total energy versus time for liquid isocyanic acid 
with a nuclear temperature of 300 K and T e = K computed 
with "exact" Born-Oppenheimer MD and the Fast QMMD 
method with exact and approximate linear scaling density 
matrix constructions. The numerical threshold r is applied in 
the linear scaling solver [26] below which all matrix elements 
are set to zero after each iteration. 



in energy is several orders of magnitude smaller than in 
previous attempts to combine linear scaling solvers with 
regular Born-Oppenheimer molecular dynamics [39T - I53] . 
The trajectories computed with an electronic temper- 
ature corresponding to k^,T e — 0.5 eV differ qualitatively 
from those computed with zero electronic temperature. 
Figure [6] and Table [n] show that while the "exact" Born- 
Oppenheimer trajectory conserves the free energy to an 
extremely high tolerance over the duration of the simula- 
tion, the total free energy in the fast quantum mechanical 
molecular dynamics simulation exhibit random-walk be- 
haviour about the mean value. Although the fast quan- 
tum mechanical molecular dynamics simulations involve 
an approximate expression for the entropy, we find that 
this alone cannot account for the level of fluctuations ob- 
served. Instead, we have found that the rescaling of the 
k value in the integration, Eq. (16), affects this random- 
walk. By changing the rescaling factor to 3/4, instead 



of 1/2 as in all the other examples, the amplitude of the 
random walk is significantly reduced. Nevertheless, the 
fast quantum mechanical molecular dynamics trajecto- 
ries at finite electronic temperature exhibit systematic 
drifts in the total energy that are negligible and the fluc- 
tuations of the total energy about the mean are of the 
same order as those that arise from the application of the 
approximate linear scaling method at T e = K. 



IV. CONVERGENCE PROPERTIES 



The fast quantum mechanical molecular dynamics 



scheme, Eqs. (14 1 and (15), can also be analyzed in 
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FIG. 6: Total free energy versus time for liquid isocyanic acid 
with a nuclear temperature of 300 K and ksT e = 0.5 eV com- 
puted with "exact" Born-Oppenheimer molecular dynamics 
(BOMD) and the fast quantum mechanical molecular dynam- 
ics (Fast-QMMD) method with k rescaled by 3/4 instead of 
1/2. 



TABLE II: Standard deviation, a, of the total energy about 
its mean value and the upper bound of the systematic drift 
of the total energy, -Edrift, computed from "exact" Born- 
Oppenheimer molecular dynamics (BOMD) and fast quantum 
mechanical molecular dynamics (Fast-QMMD) simulations of 
liquid isocyanic acid. The simulation were performed with the 
latte program in the non-orthogonal formulation, i.e. as im- 
plemented in Eqs. (|50|l, (|51[) and ( 55 1 with the entropy term 



approximated by Eq. (156 



ksT e 
(eV) 




0»v) 


-Edrift 

(/^eV/atom/ps) 


0.0 


Fast-QMMD (r = 0) 
Fast-QMMD (r = 10~ 5 ) 
BOMD (4 SCF/step) 


0.315 
0.702 
0.358 


5.10 x 10- 3 

0.285 
9.94 x 10" 3 


0.5 


Fast-QMMD (l/2)« 
Fast-QMMD (3/4)k 
BOMD (4 SCF/step) 


2.47 
0.786 
0.361 


1.43 

7.85 x 10" 2 
8.50 x 10~ 2 



FIG. 7: The root mean square deviation (RMSD) between the 
fast quantum mechanical molecular dynamics, Eqs. ( 14 1-( 15 1, 
and "exact" (4 SCF/step) Born-Oppenheimer molecular dy- 
namics, for the nuclear forces, the net Mulliken charges and 
the total energy for a Naphthalene molecule at room temper- 
ature. The simulation were performed with the latte molec- 
ular dynamics program using self-consistent-charge tight- 
binding theory in an orthogonal formulation at T e = 0, i.e. as 
implemented in Eqs. (571, (581 and (60 1 with 5 = 0.. 



der St 2 with a small pre-factor. This convergence demon- 
strates a consistency between the fast quantum mechan- 
ical scheme and Born-Oppenheimer molecular dynam- 
ics using Verlet integration, where the optimization-free 
scheme behaves as a well controlled and tunable approxi- 
mation. As in "exact" Born-Oppenheimer molecular dy- 
namics, the dominating error is governed by the local 
truncation error arising from the choice of finite integra- 
tion time step St, which is much larger than any difference 
between the fast quantum mechanical molecular dynam- 
ics and Born-Oppenheimer molecular dynamics. 



SUMMARY AND CONCLUSIONS 



terms of the convergence to "exact" Born-Oppenheimer 
molecular dynamics as a function of the finite integra- 
tion time step 5t. By comparing the deviation in forces, 
net Mulliken charges, and the total energy, between the 
fast quantum mechanical molecular dynamics scheme 
and an "exact" Born-Oppenheimer molecular dynamics 
as a function of 5t we can study the consistency between 
the two methods. Figure [7] shows the difference between 
a fully converged "exact" Born-Oppenheimer molecular 
dynamics simulation and the fast quantum mechanical 
molecular dynamics scheme as measured by the root 
mean square deviation over 200 fs of simulation time. We 
find that the deviation of the nuclear forces, the charges 
{qi}, as well as the total energy difference are of the or- 



Based on a free energy formulation of extended La- 
grangian Born-Oppenheimer molecular dynamics in the 
limit of vanishing self-consistent-field optimization, we 
have derived and demonstrated a fast quantum mechani- 
cal molecular dynamics scheme, Eqs. ( 14) and fl5|), which 



with a high precision can match the accuracy and fidelity 
of Born-Oppenheimer molecular dynamics. In addition 
to the removal of the self-consistent-field optimization we 
have also demonstrated compatibility with low pre-factor 
linear scaling solvers. The combined scheme provides a 
very efficient, energy conserving, low-complexity method 
to perform accurate quantum molecular dynamics simu- 
lations. Our findings show how the computational gap 
between classical and quantum mechanical molecular dy- 
namics simulations can be reduced significantly. 
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VII. APPENDIX 

A. Calculating the forces in Hartree-Fock theory 

Here we present some details of the fast quantum me- 
chanical molecular dynamics, Eqs. (14) and (15), using a 



simple but general Hartree-Fock formalism, which should 
be directly applicable to a broad class of hybrid and semi- 
empirical electronic structure schemes. Instead of the 
auxiliary density variable n(r) we will here use the more 
general density matrix P. In this formalism the extended 
free-energy Lagrangian in Eq. Q is given by 



£(R, R; P, P) 



\Y. M *% 



-U[TL;P}+T e S[R;P} 



(17) 



+ ^Tr[P 2 ]-^ 2 Tr[(V gs -P) 2 }, 

with the potential energy chosen as 

W[R; P] = 2Tr[hD(P)} + Tr{D(P)G[D{P)}}, (18) 

and ground state (gs) density matrix T> gs . S\R;P] is 
an unspecified electronic entropy term, which will be de- 
termined by the requirement to make the derived forces 
variationally correct, and T e is the electronic tempera- 
ture. The Fockian (or the effective single-particle Hamil- 
tonian) is 



F[P]=h + G[P], 



(19) 



with the short-hand notation, G[P] = 2J[P] - K[P], 
where J[P] and K[P] are the conventional Coulomb and 
exchange matrices, and h is the matrix of the one-electron 
part 5,6. The temperature dependent density matrix 



D(P) = z(e^ F± ^-^ + l) V, 



(20) 



which corresponds to a[n] in Eq. ([5]), is given as a Fermi 
function of the orthogonalized Fockian, 

F ± [P} = Z T F[P]Z. (21) 

Here Z and its transpose Z T are the inverse Lowdin 
or Cholesky-like factors of the overlap matrix, S, deter- 
mined by the relation 

Z T SZ = I. 



At zero electronic temperature (T e = K) the Fermi- 
operator expansion in Eq. ( 20 ) is given by the Heaviside 



step function, with the step formed at the chemical po- 
tential /iQ, separating the occupied from the unoccupied 
states. 

The Eulcr-Lagrange equations of motion of C in Eq. 



(171 are given by 



M k R k = - 



dU 
dlf k 



T e 



dS 

dRk 



2^8R k 



Tr[(V gs - P) 2 ] 



(23) 



and 



fiP = ^tu 2 (V gs -P) 



dU_ 
dP 



T, 



R 



dS_ 
dP 



(24) 



R 



A cumbersome but fairly straightforward derivation (see 
Ref. [H] for a closely related example) , using the relation 
and notation Z Rk = dZ/dR k = -(l/2)S- 1 S Rk Z, and 
defining the S term such that 



and 



T, 



T 

-L P. 



OS 



dP- 



OS 

dm 



= 2Tr 



R 



F^\D]Di. 



= 2Tr[F ± [D]DJ ik ] 



(25) 



(26) 



gives the equations of motion 

M k R k = ~2Tr[h Rk D] - Tr[DG Rk {D)] 
+Tr[{DF[D]S- 1 + S- 1 F[D]D)S R ] 



(27) 



1 



<:) 



2^dR k 



Tr[{V gs - P) 2 ] 



and 



f iP = fiuj 2 (V gs -P). 



(28) 



Notice that because of matrix symmetry Py is not in- 
dependent form Pji. The partial derivatives of matrix 
elements Py are therefore both over Py and Pji. In the 
limit fj, — > 0, we get the final equations of motion for the 
fast quantum mechanical molecular dynamics scheme, 



M k R k = -2Tr [h Rk D] - Tr [DG Rk (D)} 
+Tr [(DF^S- 1 + S- 1 F[D]D)S Rk ] , 

P = oj 2 (D(P) - P) , 



(29) 



(30) 



where we have included the substi tuti on of T> gs with 
D(P) in the same way as in Eq. (15), i.e. with us 2 



(22) rescaled by a constant c < 1. The notation for the 



partial derivative of the two-electron term is denned as 
G Rk (D) = dG(D)/dRk\ D , i.e. under the condition of 
constant density matrix D. 

The last term in Eq. ( [29] ), which includes the basis-set 
dependence Sr h is the Pulay force term that here is given 
in a generalized form that is valid also for non-idempotent 
density matrices at finite electronic tememperatures |42| . 



B. Approximate Entropy contribution 

The <S[R; P] term is defined such that the two con- 
ditions in Eqs. p5| and (26) are fulfilled. At the self- 
consistent ground state density, i.e. when P = D = T> 



gs- 



both these conditions are automatically satisfied by the 
corresponding regular ground state (gs) electronic en- 
tropy contribution to the free energy [9], 



5 gs [R;P]=«S gs [R; J D J -(P)] 

= -2k B Tr[I> L \n{D^) + (I - D^) ln(7 - D^% 



(31) 



where the relation between D 1 - and D is given by the 
congruence transformation 



D = ZD^Z T . 



(32) 



A related derivation is given in Ref. |42j. Using the 



approximate estimate <S gs [R;P] in Eq. (31) when P and 
D deviate from the ground state gives, 



T 



dS 
dS 



dR k 



= 2Tr 



R 



fHp]dI 



2Tr[F ± [P}DJ ik ]. 



(33) 



which only approximately fulfills the conditions in Eqs. 
( 25 1 and ( 26 1 . It is possible to show that the error is 
linear in SP = D — P by a linearization of F 1 - [D] around 
P. Since D(P) and P can be assumed to be close to the 
ground state, SP is small. From the scaling result illus- 
trated in Fig. [7] the error should therefore be quadratic 
in the integration time step, i.e. ~ St 2 . We may therefore 
approximate the total free energy using S gs [R; P] , which 
is zero at T e = K. However, for the exact formulation 



and derivation of the equations of motion, Eqs. (29) and 



(30), the entropy contribution, T e <S[R;P], is unknown, 
both at finite and zero temperatures. As is seen in the 
equations of motion, Eqs. (f29J) and (30), this does not 



affect the forces or the dynamics, only the estimate of 
the constant of motion, 



E 



'tot 



u 



(34) 



is approximated. By comparing the approximate E tot 
to optimized "exact" Born-Oppenheimer molecular dy- 
namics simulations, the accuracy of the dynamics can be 
estimated. 



C. Alternative potential energy forms 

As an alternative to the potential energy, U(R.;P), in 
Eq. ( 18 ) we may chose other functional forms that are 
equivalent at the ground state, i.e. when P = D = T> gs . 
By using the Harris-Foulkes-like relation [331 [M] , 

Tr[DG(D)] « Tr[(2D - P)G(P)}, 

which has an error of second order in SP = D 
may, for example, choose 

W[R; P] = 2Tr[hD{P)] + Tr{[2D{P) - P]G{P)}, (36) 

as our potential energy term. In this case, the equations 
of motion at T e = corresponding to Eqs. (29) and (30) 
become 



(35) 
P, we 



M k R k = -2Tr [h Rk D] - Tr{[2D - P}G Rk {P)} 

+Tr [(DF^S- 1 + S- 1 F[P}D)S Rk ], 
and 

P = u? (D-P), 
with the constant of motion 

E tot = \Y J M k R 2 k + 2Tr[hD] 



+Tr{{2D - P)G{P)} - T e S(R; P). 

The entropy term that makes the nuclear forces varia- 
tionally correct is here fulfilled by the expression in Eq. 
(31). With this choice of potential our dynamics only 



(37) 



(38) 



(39) 



requires one Fockian or effective single particle Hamito- 
nian construction per time step. Unfortunately, the error 
in the Pulay force has been found to be large compared 
to Eq. (p9|. The dynamics in Eqs. p7| and (1381) should 



therefore be used only for orthogonal representations, i.e. 
when the overlap matrix S = I. 



D. Self-Consistent-Charge Density Functional 
Tight-Binding Theory 

In self-consistent-charge density functional based tight- 
binding theory [4"TH50"] the continuous electronic density, 
<r(n), or the density matrix, D(P), in Eq. ( 18 ) is replaced 
by the net Mullikcn charges q[n] = {qi} for each atom i, 
where n = {n^} are the dynamical variables correspond- 
ing to P. The potential energy functional U in Eq. ( 18 ) 
is then reduced to 



W[R;n] = 2 ^ e i --^ gi (n)g j (n)^+S patr [R]. 



(40) 



Here Si are the (doubly) occupied eigenvalues of the 
charge dependent effective single-particle Hamiltonian 

(41) 



-(!/ 2 ) 2^ (Sia,kf3>V k e p, Jl3 + V^ kj3 
k/3> 



Vf^ ha ,S k /3' y j/3) 
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where 



^"*0' = S«( n )7iJW, 



P'Pi 



(42) 



hiajp is a parameterized Slater-Koster tight-binding 
Hamiltonian, Si a jp the overlap matrix, i and j are 
atomic indices and a and f3 are orbital labels [51] , The 
net Mulliken charges are given by 



qi[n] = 2^ (el atia - Q° iaA 



(43) 



with the density matrix 

g 1 - = g ± [n] = (g^^H-/*/) + ^ _1 _ (44) 

using the orthogonalized Hamiltonian 



ff-^n] = Z T iJ[n]Z. 



(45) 



Here g° is the density matrix of the corresponding sepa- 
rate non-interacting atoms. The de-orthogonalized den- 
sity matrix is 



g = g[n] = Zq- l [vl]Z' 1 



(46) 



and as above, the congruence transformation factors are 
defined through 



Z T SZ = /, 
where S is the basis set overlap matrix. 



(47) 



The electron-electron interaction in Eq. (40) is deter- 



mined by 7jj, which decays like 1/R at large distances 
and equals the Hubbard repulsion for the on-site interac- 
tion. _E pa i r [R] is a sum of pair potentials, <f>(R), that pro- 
vide short-range repulsion. The radial dependence, £(R), 
of the Slater-Koster bond integrals, elements of the over- 
lap matrix, and the 4>{R) are all represented analytically 
in LATTE by the mathematically convenient form, 



C(R)=Aol[exp(A i R i 



(48) 



i=l 



where A$ to A4 are adjustable parameters that are fitted 
to the results of quantum chemical calculations on small 
molecules. To ensure that the off-diagonal elements of h 
and S and the 4>(R) in our self-consistent tight-binding 
implementation decay smoothly to zero at a specified dis- 
tance, i?cut, we replace the C(R) by cut-off tails of the 
form, 



t(R) = B Q + AR(Bi + AR(B 2 

+ AR(B 3 + AR{B A 



ARB 5 )))) (49) 



at R = R\, where AR = R — i? x and B to B 5 are 
adjustable parameters. The adjustable parameters are 
parameterized to match the value and first and second 
derivatives of t(R) and ((R) at R = R\ and to set the 
value and first and second derivatives of t(R) to zero at 

R = Rent- 



1. Non-orthogonal representation at T e > 

The fast quantum mechanical molecular dynamics 
scheme, Eqs. J29} and Q or Eqs. (JsJ) and ^, using 
self-consistent tight-binding theory in its non-orthogonal 
formulation is given by 



M k R k = -2Tr[gH Rk ] 

1 >-^ d^n >-^ dqj 



*>j 



1 -j 



(50) 



+Tr[(S- 1 H[q]g+gH[q]S- 1 )S Rk 
and 



dE pail [R] 
dR k ' 



hi = uj 2 (q t - n t ) 



where 



and 



Hr^ = 



St 



dH 

dRk 



OS 
dRk 



(51) 



(52) 



(53) 



The partial derivatives of qj and H in Eqs. (50) and 



( 52 ) are with respect to a constant density matrix g in 



its non-orthogonal form, i.e. including an S dependence 
ofq-j, 



dqj 



= 2 22(eS Rk ) jada . 

3 a£j 



dR k 
The total energy is given by 

Etot = lJ2 M * k k+2j2 £ * 



(54) 



ifEocc 



ij2 WW + £ P air[R] - T e S[R; n], 



(55) 



1 -j 



with the entropy contribution to the free energy approx- 
imated by 

5[R; n] « -2fc B J2 if* HA) + ^ ~ & h < 1 ~ /*)> ■ 

(56) 
Here /j = ft [n] are the eigenstates of the Fermi operator 
expansion g 1 - [n] of H 1 - [n] in Eq. ( 44 ) . 



2. Orthogonal representation at T e — 

For orthogonal formulations, i.e. when S = I, and at 
zero electronic temperature, T e = 0, we will base our 



11 



dynamics on the equations of motion in Eqs. (37) and 



38| . In this case the fast quantum mechanical molecular 



dynamics scheme, Eqs. ( 50 )-( 51 ), is given by 



M k R k = -2Tr[gH Rk ] 



<9£pair [R] 

dR k 



2^ 

i,3 



ran. 



dR k 



hi = u 2 (qi - m) , 



(57) 



(58) 



where 



{H Rk [xi]} ia jp = 






0R k 



Otik 



The density matrix is given directly from the step func- 
tion of the Hamiltonian, g — 0(/iqI — H [n]), without any 
de-orthogonalization that requires the calculation of the 
inverse factorization of the overlap matrix, Eq. ( 47 1 . The 
constant of motion, E to ti is approximate by 



E tot ^ 1 -Y,M k R 2 k + 2j2^ 

k iGocc 

--^(2tii- Qi)njjij + E pail [R]. 



(60) 



3. General remarks 



Apart from the first few initial molecular dynam- 
ics time steps, where we apply a high degree of self- 
consistent-field convergence and set n — q, no ground 
state self-consistent-field optimization is required. The 
density matrix, q, and the Hamiltonian, H, necessary in 
the force calculations (and for the total energy) are calcu- 
lated only once per time step in the orthogonal case with 
one additional construction of the Hamiltonian required 
in non-orthogonal simulations. The numerical integra- 
tion of the equations of motion in Eq. ( 14 ) is performed 



with the velocity Verlet scheme and in Eq. (151 with the 



modified Verlet scheme in Eq. ( 16 ) as described in Ref. 
|36) . For the examples presented here we used the mod- 
ified Verlet scheme including dissipation (a > 0) with 
K = 5 and the constant k — St 2 oj 2 as given in Ref. [31)] 
was rescaled by a factor 1/2 in all examples except for 
one of the test cases in Fig. [6| 
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